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Abstract 

We study a model for a two-mode atomic-molecular Bose-Einstein condensate. 
Starting with a classical analysis we determine the phase space fixed points of the 
system. It is found that bifurcations of the fixed points naturally separate the 
coupling parameter space into four regions. The different regions give rise to quali- 
tatively different dynamics. We then show that this classification holds true for the 
quantum dynamics. 
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1 Introduction 



Since the early experimental realisations of Bose-Einstein condensates (BECs) using al- 
kali atoms [1,2], a significant effort has been made to produce a stable BEC in a gas of 
molecules [3]. A molecular condensate could lead to a host of new scientific investiga- 
tions that includes the quantum gas with anisotropic dipolar interactions [4], the study 
of rotational and vibrational energy transfer processes [5] and coherent chemistry where 
the reactants and products are in a coherent quantum superposition of states [6] , among 
others. In recent years the creation of a molecular BEC from an atomic BEC has been 
achieved by different techniques such as photoassociation [7], two-photon Raman transi- 
tion [8] and Feshbach resonance [9]. 

Prom a theoretical point of view, molecular BECs may be studied using the Gross- 
Pitaevski (GP) equations and mean-field theory (MPT) (e.g. see [10,11]). The GP-MPT 
approach reduces the full multi-body problem into a set of coupled nonlinear Schrodinger 
equations, which are then solved numerically to obtain the Josephson-type dynamics 
of the coupled atomic and molecular fields. An approximation can be made to reduce 
the complex multi-body problem into a two-mode problem. An analysis of this two- 
mode Hamiltonian was carried out in [10], where it was established that the quantum 
solutions break away from the MPT predictions in the vicinity of the dynamically unstable 
molecular mode due to strong quantum fiuctuations. It has been shown that the two- 
mode Hamiltonian is an exactly solvable model in the framework of the algebraic Bethe 
ansatz method [12] and an analysis using these results was given in [13]. However in most 
of the above investigations, the atom-atom, atom-molecule and molecule-molecule 5'-wave 
scattering interactions were not taken into account. 

In the present work we focus on a more general Hamiltonian which takes into account 
the S'-wave scattering interactions. By means of a classical analysis we first obtain the 
fixed points of the system and find that the space of coupling parameters divides into four 
distinct regions which are determined by fixed point bifurcations. By contrast, only three 
such regions exist when the S-wscve scattering interactions are neglected. The results allow 
us to qualitatively predict the dynamical behaviour of the system in terms of whether the 
evolution is localised or delocalised. Using exact diagonalisation of the Hamiltonian, we 
then see that the quantum dynamics within each region has a similar character. 

The paper is organised as follows: In section 2 we present the Hamiltonian and in 
section 3 a classical analysis of the model is performed. In section 4 we investigate the 
quantum dynamics through the time evolution of the expectation value of the relative 
atom number. Section 5 is reserved for a discussion of the results. 



2 The model 

Let us consider the following general Hamiltonian, based on the two-mode approximation, 
describing the coupling between atomic and diatomic-molecular Bose-Einstein conden- 
sates 

H = UaNl + U^Nl + UabNaN, + HaNa + l^bN, + n{o)o)h + h'^ aa) . (1) 
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Above, is the creation operator for an atomic mode while 6"^ creates a molecular mode. 
The Hamiltonian commutes with the total atom number N = Na + 2Nh, where Na = a"^a 
and Nh = b'^b. Notice that the change of variable Q — ^7 is equivalent to the unitary 
transformation 

b -b. (2) 

The parameters Uj describe S-wave scattering, taking into account the atom-atom 
(Ua), atom- molecule (Uab) and molecule- molecule {Uf,) interactions. The parameters 
are external potentials and Q is the amplitude for interconversion of atoms and molecules. 
In the limit Ua = Uab = = 0, (UJ has been studied using a variety of methods [10,12-14]. 
However in the experimental context, the 5*- wave scattering interactions play a significant 
role. It will be seen below that for the general model the inclusion of these scattering 
terms has a non-trivial consequence. We mention that generally the values for Ub and Uab 
are unknown [6,8], although some estimates exist in the case of ^^Rb [11]. 

We finally note that the Hamiltonian (P) is a natural generalisation of the two-site 
Bose-Hubbard model 

H = U{N, - N^f + ^l{Nl - N2) + n{a\a2 + ala,) (3) 

which has been extensively studied as a model for quantum tunneling between two single- 
mode Bose-Einstein condensates [14-20]. Our analysis will show that despite apparent 
similarities between the Hamiltonians (^Q) and 0, they do display some very different 
properties. This aspect will be discussed in Sectional 

3 The classical analysis 

Let Nj, 6j, j = a, b be quantum variables satisfying the canonical relations 

%, Ob] = [Na, Nb] = 0, [N„ Ok] = 

Using the fact that 

exp(i6'j)A/'j = {Nj + l)exp(26'j) 
we make a change of variables from the operators j, j\ j = a, b via 

j = exp{iej)^/N'j, = a/A^ expi-iOj) 

such that the Heisenberg canonical commutation relations are preserved. We make a 
further change of variables 

z = j^iNa-2Nb), 

N = Na + 2Nb, 
N 

= ^(2^a — Ob), 

such that z and 6 are canonically conjugate variables; i.e. 

[z, 6] = il. 
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In the limit of large we can now approximate the (rescaled) Hamiltonian by 

H = \z' + 2az + (3 + ^2{l ~z){l + z) cos (4) 



with 
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where, since is conserved, we treat it as a constant. We note that the unitary transfor- 
mation (j2)) is equivalent to ^ O + Ntc/A. Also, since the Hamiltonian is time- reversal 
invariant, we will hereafter restrict our analysis to the case A > 0. 

We now regard (jU as a classical Hamiltonian and investigate the fixed points of the 
system. The first step is to find Hamilton's equations of motion which yields 

dz OH 4 /— , . [AO 

dO OH ^ l-Sz [Ad 

— r = -TT- = 2Xz + 2a H , cos — 

dt dz ^2(1 -2) 

The fixed points of the system are determined by the condition 

dH dH 

Due to periodicity of the solutions, below we restrict to ^ G [0, Nti/2). This leads to the 
following classification: 

• 9 = Ntt/A, and z is a solution of 

l-3z 

Xz + a 



2^2(1-^) 

which has no solution for X — a < —1 while there is a unique locally minimal solution 
for A — a > —1. 

= 0, and 2; is a solution of 

A 3z — 1 

Xz + a = , (o) 

2^/2(1^ 

which has a unique locally maximal solution for A — a < 1 while for A — a > 1 
there are either two solutions (one locally maximal point and one saddle point) or 
no solutions. In Fig. [T]we present a graphical solution of (jH)). 
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• z = —1 and 6' is a solution of 




= X — a 



for which there are two saddle point solutions for |A — q;| < 1. 

It is also useful to identify the points z = 1, 6 = Ntt/8 and z = 1, 6 = SNn/S, where the 
singular derivative dH/dz changes sign. For z = 1 the Hamiltonian ^ is independent 
of 9, so these points essentially behave like a saddle point. We remark that (jlj) is also 
independent of 6 for z = —1. 




z z 



Figure 1: Graphical solution of equation The crossing between the straight line (left 
hand side of eq.®) and the curve (right hand side of eq.®) for different X — a values 
represents the solution(s) for each case. There is just one solution on the left (A — a < 1), 
while there are either two solutions or no solution on the right (A — a > 1). 

From the above we see that there exist fixed point bifurcations for certain choices of 
the coupling parameters. These bifurcations allow us to divide the parameter space into 
four regions, as depicted in Fig. |21 The asymptotic form of the boundary between regions 
I and II is discussed in the Appendix. 

To visualise the dynamics, it is useful to plot the level curves of the Hamiltonian 
dH). Since the fixed point bifurcations change the topology of the level curves, qualitative 
differences can be observed between each of the four regions. The results are shown in 
Figs. (|3l4p ■ where for clarity we now take iO/N E [—2%, 27i]. 

Fig. Efa) shows the typical character of the level curves in region I. The maximal level 
curve occurs along the phase space boundary z = —1 and there are two local minima. 
Note that for no choice of parameters do these minima occur on the boundary z = 1, but 
they may occur arbitrarily close to this boundary. If the initial state of the system has 
z ^ 1 then z will remain close to 1 for all subsequent times. A similar situation is true 
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Figure 2: Parameter space diagram identifying the different types of solution for equation 
(0). In region I there are no solutions for z when ^ = 0, and one solution for z when 
Q = Nn/A. In region II there are two solutions for z when ^ = 0, and one solution for z 
when 6 = Nn/A. In region III exists one solution for z when 6 = 0, one solution for z 
when 6 = Nn/A, and two solutions for 6 when z = —1. In region IV there is one solution 
for z when ^ = 0, and no solution for z when 6 = Ntt/A. The boundary separating regions 
II and III is given by A = a + 1, while the equation X = a — 1 separates the regions III 
and IV. The boundary between regions I and II has been obtained numerically. 



if the initial state of the system has z ~ — 1. For both cases we see that the evolution of 
the system is localised. 

As the coupling parameters are changed and the system crosses the boundary into 
region II, two new fixed points, a maxima and a saddle point, emerge at 6' = which can 
happen for any z G [—1, 1). On crossing this parameter space boundary the maximum 
may move towards the phase space boundary z = 1 while the saddle point approaches 
z = —1. Also the two minima may move away from the phase space boundary z = —1, as 
depicted in Fig. ISfb). The consequence for the dynamics is that for an initial state with 
z ^ —1 the evolution of the system is still localised, but for an initial state with z ^ 1 
the evolution is delocalised. 

Fig. Ilfa) illustrates what happens when the coupling parameters are tuned to cross 
over from region II into region III. The saddle point at 6 = approaches z = —1, reaching 
the phase space boundary exactly when the coupling parameters lie in the boundary 
between regions II and III. The saddle point then undergoes a bifurcation into two saddle 
points occurring at ^ = —1 for different values of 6 in region III. The two mimima have 
also moved away from z = 1 towards z = —1. Now the dynamics is delocalised for both 
initial states z ^ 1 and z ^ —1. 

It is also possible to tune the parameters to move from region I directly to region III. 
At the boundary between the regions, a single local maximum emerges from the point 
z = —1, 6 = 0. As the parameters are tuned to move away from the boundary into region 
III, the maximum moves towards z = 1 while the minima at ^ = ±iV7r/4 approaches 
z = -1. 

Moving from region III towards region IV causes the two saddle points for z = —1 to 
move towards 6 = ±A^7r/4. Again, the two minima for 6 = ±Ntt/4 move towards z = —1. 
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Figure 3: Level curves of the Hamiltonian in (a) region I and (b) region II. The 
parameter values are A = 1.0, a = —8.0 for region I and A = 1.0, a = —0.2 for region II. 
In region I we observe the presence of local minima for W/N = ivr. Besides the minima 
at 46/N = ±7r, two additional fixed points (a maximum and a saddle point) are apparent 
in region II occurring at ^ = 0. 

Each minimum converges with a saddle point exactly when the coupling parameters are 
on the boundary of regions III and IV. Varying the coupling parameters further into 
region IV we find that minima for the Hamiltonian are always at z = —1, 6 = ±A^7r/4, 
and the local maximum for ^ = lies close to 2; = 1, as shown in Fig. For this case the 
dynamics is localised for both initial states z ^ 1 and z ^ —1. 

The above discussion gives a general qualitative description of the dynamical behaviour 
of the classical system in terms of the four regions identified in the parameter space. We 
emphasise that the change in the classical dynamics as the boundary between two regions 
is crossed is smooth. Nonetheless, the analysis does give a useful insight into the possible 
general dynamical behaviours. Below we will show that the same holds true for the 
quantum dynamics. 

4 Quantum dynamics 

Having analysed the classical dynamics, we now want to investigate the extent to which 
a similar scenario holds for the quantum system. For the case A = (where the coupling 
for all S-wave scattering interactions is zero) the quantum dynamics has previously been 
studied in [10,13]. In this instance region II is not accessible. It was shown that the 
dynamics is delocalised for |a| < 1 and localised otherwise for both atomic and molecular 
inital states, consistent with the classical results described above. A surprising aspect of 
the classical analysis is the existence of region II where the evolution of a purely molecular 
inital state is highly locahsed, whereas the evolution of a purely atomic initial state is 
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Figure 4: Level curves of the Hamiltonian (0} in (a) region III and (b) region IV. The 
parameter values are A = 1.0, a = 0.2 on the left and A = 1.0, a = 3.0 on the right. In 
region III we observe the presence of minima at 46'/A^ = ±7r and for ^ = just one fixed 
point, a maximum. There are also saddle points for when z = —1. In region IV just one 
fixed point (a maximum) occurs for ^ = 0, which always has z < 1. In contrast the global 
minimum occurs for z = —1. 

completely delocalised. We will see that this also occurs for the quantum case. Thus the 
inclusion of the S'-wave scattering interactions into the Hamiltonian gives richer dynamics. 

The time evolution of any state is given by = U{t)\(j)), where U{t) is the temporal 

operator U(t) = ^^^q |m) (m| exp(— ii^m^), 1"^-) is an eigenstate with energy Em and 
10) = \Na, Nh) represents the initial Fock state with A^^ atoms and A^^^ molecules such that 
Na + 2Ni, = N. We adopt the method of directly diagonalising the Hamiltonian as done 
in [19, 20] for the Bose-Hubbard Hamiltonian Q and compute the expectation value of 
the relative number of atoms 

(Nait) - 2N,{t)) = {^{t)\N, - 2N,\^{t)) 

using two different initial state configurations: a purely atomic state and a purely molecu- 
lar state. Hereafter, we will fix the following parameters N = 100, Q = 1.0, fia = 0.0, fib = 
0.0 and Ub = 1.0. 

In Fig. |31we plot the expectation value of the relative number of atoms for A = 1.0 
and the choices a = -8.0, -0.2, 0.2, 3. The graphs depict the quantum dynamics for those 
cases where the system is in regions I, II, HI and IV from top to bottom respectively. 
On the left we are using a purely atomic initial state |A^, 0) and on the right hand side a 
purely molecular initial state |0, A^/2). 

Figure El displays aspects of the quantum dynamics, such as the collapse and revival 
of oscillations and non-periodic oscillations, which are not features of the corresponding 
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Figure 5: Time evolution of the expectation value of the imbalance population {Na — 
2Nf,)/N in the four regions defined by the diagram with a purely atomic initial state 
|A^, 0) on the left and a purely molecular initial state |0, on the right. We are using 
A = 1.0 and a = —8.0, —0.2, 0.2, 3.0 (or, in terms of the original variables, Ub = l,Ua = 
-0.881, 0.222, 0.278, 0.674 and Uab = -1.546, 0.660, 0.774, 1.566). 



classical dynamics (cf. [15] for analogous results for the Hamiltonian (jHl)). However it 
also shows that the classification based on classical fixed point bifurcations to determine 
whether the dynamic evolution is localised or delocalised applies to the quantum case. In 
particular, in region II it is clear that for an initial atomic state the evolution is completely 
delocalised, but localised for an initial molecular state. 



5 Discussion 

Using the classical Hamiltonian Q as an approximation to the quantum Hamiltonian 
(Q, we have undertaken an analysis to determine the fixed points of the system. The 
bifurcations of the fixed points divide the coupling parameter space into different regions 
characterising different dynamics, which can also be seen for the quantum dynamics. It 
is necessary to establish the extent to which the classical approximation is valid. Since A 
and a vary with the number of particles, it is required that the gap between successive 
energy levels should approach a continuum for large A^. This imposes that A, a « N^/^. 
We can compare this situation to the case of the Bose-Hubbard model 0, where a similar 
classical analysis is valid for |f//^^| « N [17]. It was shown in [20] that for that model 
there are transitions in the dynamical behaviour for the quantum regime >> A^, 

which are not apparent from the classical analysis. These properties were found to be 
closely related to couplings for when the energy gap between the ground and first excited 
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state was minimal or maximal. We should expect a similar result to occur for (^. 

The relationship between fixed point bifurcations and ground-state entanglement has 
been studied in [25]. There it was argued that the ground state entanglement of a quantum 
system will be maximal whenever the classical system undergoes a supercritical pitchfork 
bifurcation for the lowest energy in phase space. A peak in a measure of ground-state 
entanglement has been shown in many instances to be indicative of a quantum phase 
transition [21-23]. For the Hamiltonian (P) we have considered here, there are no super- 
critical pitchfork bifurcations. For A = there is a quantum phase transition at a = 1, as 
can be seen from the behaviour of certain ground state correlation functions [11,13,14]. 
This does correspond to a bifurcation of the lowest energy in phase space. Calculation of 
the ground-state entanglement in this case have been undertaken in [14] showing that it 
is maximal at a coupling different from the critical point. 

This is in some contrast to the Bose-Hubbard model Q. There, a supercritical pitch- 
fork bifurcation of the lowest energy occurs in the attractive case [18,20], and the results 
of [24] suggest that indeed the entanglement is maximal at this coupling. (For the repulsive 
case the ground state entanglement is a smooth, monotonic function of the coupling [14].) 
However the transition from localisation to delocalisation for the dynamics as studied 
in [15,19,20] does not occur at the bifurcation. Despite the apparent similarities between 
(0) and Q, we can see that the inter-relationship between bifurcations of the classical 
system and properties of the quantum system are very different. 
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Appendix 

In this Appendix we analyse the boundary dividing regions I and II. In particular, we 
determine the asymptotic relation between A and a when A is large and when A is close 
to 1/2. We also compute the maximum value of a on this boundary. 



where the fixed points occur when f{z) = g{z). We want to determine the boundary 
between the cases when there is no solution and two solutions. This boundary is given 
by the case when f{z) is the tangent line to g{z). Now 
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^ = ^(3(1 - ^ _ _ )-3/2) 

so z is determined by the condition 

Below we consider three cases: 

(i) First put z = — 1 + where u is small and positive. Then 

dg 1 3 

— ~ — \ u. 

dz 2 16 



Solving for u gives 



Now we need 



_(2A- 1). 



X{-l + u) + a 
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We can substitute in u to find a relation between A and a: 

1 A 1\ 16 1 " 
Q,^ ^ I X I X 

2V 27 3 V 2 
This curve is vahd for (A — 1/2) positive. Also 

= 1 



da 
dX 



A=l/2 



SO the curve separating regions II and III is tangential to the curve separating regions 
I and III at A = 1/2. 

(ii) Next we look at the case when z — 1 — u with u small and positive. Here we find 



2-3u 



2V2u 
1 



2u 
1 



dz ^ 2v^ti3/2 
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so that 

2 

This leads to 



a = g{z) — Xz 

A + (8) 

The asymptotic equation (jSJ is vahd for large positive values of A. 

(iii) To complete the picture, finally we investigate the maximum of a with respect to 
A. From ()7|8p we have 

da dg dz ^dz 

dX dz dX dX 

= -z (10) 

so the maximum occurs at z = 0. Looking at the asymptotic behaviour around 
2; = we have 



which gives 



qlz) ~ ^ 1 z z' 

dg 5 9 

— ~ ^ IH z) 

dz 4^2 10 ' 



10 4^2 

z ~ — A — L 

9^5 



2^ 



Using this we can find an expression for a in terms of A: 

1 25 ,4^2^ ,2 

a ~ -= -={ A - 1)^ 

2v^ 36^2 5 

The first term above corresponds to the maximal value of a ~ —0.35 as depicted in 
Fig. El 
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